# Imports
import pandas as pd
import numpy as np
from pandas.api.types import CategoricalDtype
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.graph_objects as go
from scipy import stats
import statsmodels.api as sm
from scipy.stats import bootstrap
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
from tqdm import tqdm
import kagglehub
import os

# Aesthetic configuration for plots
sns.set_theme(style="whitegrid")
plt.rcParams['figure.figsize'] = (10, 6)

# Download the latest version of the dataset
path = kagglehub.dataset_download("pratyushpuri/global-ai-job-market-trend-2025")

print("Path to dataset files:", path)

# Exact name of the CSV file
csv_file_name = "ai_job_dataset.csv"

# Combine the folder path with the CSV file name
full_file_path = os.path.join(path, csv_file_name)

# Ensure that the CSV file has been downloaded and loaded correctly
try:
    # Read the CSV file using the full path
    # Apply 'parse_dates' to the date columns
    df = pd.read_csv(
        full_file_path,
        parse_dates=['posting_date', 'application_deadline']
    )
    print(
        f"✅ Dataset '{csv_file_name}' successfully loaded: "
        f"{df.shape[0]} rows, {df.shape[1]} columns."
    )
except FileNotFoundError:
    print(
        f"⚠️ ERROR: '{csv_file_name}' not found. "
        "Please make sure the download process completed correctly."
    )

Path to dataset files: C:\Users\tayfu\.cache\kagglehub\datasets\pratyushpuri\global-ai-job-market-trend-2025\versions\1
✅ Dataset 'ai_job_dataset.csv' successfully loaded: 15000 rows, 19 columns.

import plotly.io as pio
# This option forces Plotly to be saved as a full HTML block inside the notebook file
pio.renderers.default = "notebook"

# Quick overview of the DataFrame
df.head()

f,c = df.shape
print(f"The DataFrame contains {f} rows and {c} columns")

The DataFrame contains 15000 rows and 19 columns

df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 15000 entries, 0 to 14999
Data columns (total 19 columns):
 #   Column                  Non-Null Count  Dtype         
---  ------                  --------------  -----         
 0   job_id                  15000 non-null  object        
 1   job_title               15000 non-null  object        
 2   salary_usd              15000 non-null  int64         
 3   salary_currency         15000 non-null  object        
 4   experience_level        15000 non-null  object        
 5   employment_type         15000 non-null  object        
 6   company_location        15000 non-null  object        
 7   company_size            15000 non-null  object        
 8   employee_residence      15000 non-null  object        
 9   remote_ratio            15000 non-null  int64         
 10  required_skills         15000 non-null  object        
 11  education_required      15000 non-null  object        
 12  years_experience        15000 non-null  int64         
 13  industry                15000 non-null  object        
 14  posting_date            15000 non-null  datetime64[ns]
 15  application_deadline    15000 non-null  datetime64[ns]
 16  job_description_length  15000 non-null  int64         
 17  benefits_score          15000 non-null  float64       
 18  company_name            15000 non-null  object        
dtypes: datetime64[ns](2), float64(1), int64(4), object(12)
memory usage: 2.2+ MB

def info_df(df):
    return pd.DataFrame({
        'Column': df.columns,
        'Non-Null': df.notnull().sum().values,
        'Null': df.isnull().sum().values,
        'Python Type': df.dtypes.values,
        'Num. Unique Values': [len(df[col].unique()) for col in df.columns],
        'Sample Values': [df[col].unique()[0:(min(5, len(df[col].unique())))] for col in df.columns],
    })

info = info_df(df)
info

df.describe(percentiles=[0.25, 0.5, 0.75, 0.95, 0.99]).T

df.describe(include="O").T

# Convert NaN to an empty string (we don't have any, but we include this as a best practice)
skills_series = df['required_skills'].fillna('')

# Split by commas -> converts each row into a list
skills_lists = skills_series.str.split(',')

# Strip whitespace around each skill
skills_lists = skills_lists.apply(lambda lst: [s.strip() for s in lst if s.strip() != ''])

# Explode -> each skill becomes its own row
skills_exploded = skills_lists.explode()

# Count frequencies
top_skills = skills_exploded.value_counts().head(10)

print("Top most in-demand skills:")
print(top_skills)

Top most in-demand skills:
required_skills
Python        4450
SQL           3407
TensorFlow    3022
Kubernetes    3009
Scala         2794
PyTorch       2777
Linux         2705
Git           2631
Java          2578
GCP           2442
Name: count, dtype: int64

from pandas.api.types import CategoricalDtype

# Define the logical ordering for each ordinal variable
size_order = CategoricalDtype(categories=['S', 'M', 'L'], ordered=True)

# Assume this order: Entry -> Mid -> Senior -> Executive
exp_order = CategoricalDtype(categories=['EN', 'MI', 'SE', 'EX'], ordered=True)

# Assume this academic order
edu_order = CategoricalDtype(categories=['Associate', 'Bachelor', 'Master', 'PhD'], ordered=True)

# Apply the transformations
df['company_size'] = df['company_size'].astype(size_order)
df['experience_level'] = df['experience_level'].astype(exp_order)
df['education_required'] = df['education_required'].astype(edu_order)

nominal_cols = [
    'employment_type', 'salary_currency', 'industry',
    'company_location', 'employee_residence', 'company_name'
]

for col in nominal_cols:
    df[col] = df[col].astype('category')
    # Since we do not specify 'ordered=True', Pandas treats these as unordered categories

info = info_df(df)
info

for col in df.columns:
    if len(df[col].unique())<45:
        print(df[col].value_counts())
        print("="*40)

job_title
Machine Learning Researcher    808
AI Software Engineer           784
Autonomous Systems Engineer    777
Machine Learning Engineer      772
AI Architect                   771
Head of AI                     765
NLP Engineer                   762
Robotics Engineer              759
Data Analyst                   759
AI Research Scientist          756
Data Engineer                  749
AI Product Manager             743
Research Scientist             742
Principal Data Scientist       734
AI Specialist                  728
ML Ops Engineer                725
Computer Vision Engineer       724
Data Scientist                 720
Deep Learning Engineer         718
AI Consultant                  704
Name: count, dtype: int64
========================================
salary_currency
USD    11957
EUR     2314
GBP      729
Name: count, dtype: int64
========================================
experience_level
MI    3781
EX    3760
SE    3741
EN    3718
Name: count, dtype: int64
========================================
employment_type
FT    3812
FL    3758
CT    3721
PT    3709
Name: count, dtype: int64
========================================
company_location
Germany           814
Denmark           778
Canada            769
France            769
Austria           765
Singapore         764
China             763
India             754
Sweden            752
Israel            751
Ireland           750
Switzerland       746
Finland           733
Japan             733
Australia         732
Netherlands       731
United Kingdom    729
United States     724
South Korea       722
Norway            721
Name: count, dtype: int64
========================================
company_size
S    5007
L    4998
M    4995
Name: count, dtype: int64
========================================
employee_residence
Sweden            790
France            781
Denmark           777
Austria           776
India             772
Germany           769
South Korea       763
Canada            762
China             761
Netherlands       758
United Kingdom    750
Switzerland       748
Ireland           740
Singapore         740
Israel            731
Australia         730
Norway            726
United States     716
Finland           710
Japan             700
Name: count, dtype: int64
========================================
remote_ratio
0      5075
50     5005
100    4920
Name: count, dtype: int64
========================================
education_required
Bachelor     3789
Associate    3785
Master       3748
PhD          3678
Name: count, dtype: int64
========================================
years_experience
0     1890
1     1828
4     1295
3     1247
2     1239
7      769
5      757
6      753
9      742
8      720
16     403
15     392
13     391
10     384
19     378
11     373
14     364
17     363
12     362
18     350
Name: count, dtype: int64
========================================
industry
Retail                1063
Media                 1045
Automotive            1020
Consulting            1020
Technology            1011
Real Estate           1007
Government             998
Healthcare             997
Telecommunications     997
Transportation         997
Finance                984
Energy                 976
Gaming                 967
Manufacturing          962
Education              956
Name: count, dtype: int64
========================================
company_name
TechCorp Inc                  980
Cognitive Computing           972
AI Innovations                964
Digital Transformation LLC    961
Future Systems                960
Quantum Computing Inc         960
Cloud AI Solutions            951
Predictive Systems            947
Smart Analytics               927
Advanced Robotics             925
Machine Intelligence Group    922
Neural Networks Co            922
Autonomous Tech               918
DataVision Ltd                909
DeepTech Ventures             897
Algorithmic Solutions         885
Name: count, dtype: int64
========================================

df.duplicated().sum()   # no duplicates

0

# Show only the first (most frequent) value in each column, along with its count
header = f"{'Column':<25} {'Mode':<30} {'Frequency':<15} {'% of non-nulls':<15}"
print(header)
print("-" * len(header))

for col in df.columns:
    column_non_null = df[col].dropna()  # drop nulls before counting rows (we don't have any, but it's good practice)
    frequencies = df[col].value_counts()
    mode_value = frequencies.index[0] if not df[col].mode().empty else None
    non_null_count = df[col].notnull().sum()
    mode_count = frequencies.values[0]
    mode_percentage = (mode_count / non_null_count * 100) if non_null_count > 0 else 0
    print(f"{col:<25} {str(mode_value):<30} {mode_count:<15} {round(mode_percentage, 2):<15}")

Column                    Mode                           Frequency       % of non-nulls 
----------------------------------------------------------------------------------------
job_id                    AI00001                        1               0.01           
job_title                 Machine Learning Researcher    808             5.39           
salary_usd                67253                          4               0.03           
salary_currency           USD                            11957           79.71          
experience_level          MI                             3781            25.21          
employment_type           FT                             3812            25.41          
company_location          Germany                        814             5.43           
company_size              S                              5007            33.38          
employee_residence        Sweden                         790             5.27           
remote_ratio              0                              5075            33.83          
required_skills           Python, TensorFlow, PyTorch    17              0.11           
education_required        Bachelor                       3789            25.26          
years_experience          0                              1890            12.6           
industry                  Retail                         1063            7.09           
posting_date              2024-07-05 00:00:00            51              0.34           
application_deadline      2025-01-05 00:00:00            47              0.31           
job_description_length    1519                           19              0.13           
benefits_score            9.9                            338             2.25           
company_name              TechCorp Inc                   980             6.53

# Code to get numeric columns
numeric_columns = df.select_dtypes(include=['number']).columns.tolist()

header = f"{'Column':<30} {'Mean':<10} {'Median':<10}"
print(header)
print("-" * len(header))
for col in numeric_columns:
    mean = df[col].mean()
    median = df[col].median()
    print(f"{col:<30} {round(mean,2):<10} {round(median,2):<10}")

Column                         Mean       Median    
----------------------------------------------------
salary_usd                     115348.97  99705.0   
remote_ratio                   49.48      50.0      
years_experience               6.25       5.0       
job_description_length         1503.31    1512.0    
benefits_score                 7.5        7.5

# Automatic alternative using bins='auto'
# 'auto' usually uses Freedman–Diaconis when the dataset is large (we have ~15,000 rows)
sns.histplot(df['salary_usd'], bins='auto', kde=True)
plt.title('Salary distribution')
plt.xlabel('Salary (USD)')
plt.show()

# Manual alternative using Freedman–Diaconis
# Compute the interquartile range (IQR)
q1 = df['salary_usd'].quantile(0.25)
q3 = df['salary_usd'].quantile(0.75)
iqr = q3 - q1

# Compute bin width (Freedman–Diaconis rule)
# h = 2 * IQR * n^(-1/3)
data_len = len(df['salary_usd'])
bin_width = 2 * iqr * (data_len ** (-1/3))

# Compute the number of bins
data_range = df['salary_usd'].max() - df['salary_usd'].min()
num_bins = int(data_range / bin_width)

print(f"Optimal number of bins (Freedman–Diaconis): {num_bins}")

# Plot using the computed number of bins
sns.histplot(df['salary_usd'], bins=num_bins, kde=True)

plt.title('Salary distribution')
plt.xlabel('Salary (USD)')
plt.show()

Optimal number of bins (Freedman–Diaconis): 59

plt.figure(figsize=(8, 6))
sns.boxplot(x=df['salary_usd'], color='lightgray')
plt.title('Salary distribution (USD)', fontsize=14)
plt.xlabel('Salary (USD)')
plt.show()

# --- ANALYSIS OF PREDICTOR VARIABLE: experience_level ---

plt.figure(figsize=(10, 5))

# Univariate plot
el_count = df['experience_level'].value_counts()

sns.barplot(x=el_count.index, y=el_count.values)
plt.title('Frequency by experience level')
plt.xlabel('Experience level')
plt.ylabel('Employees')

ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

for i, v in enumerate(el_count.values):
    plt.text(i, v + 10, str(v), ha='center', va='bottom')

plt.show()

# Bivariate plot (X vs Y)
plt.figure(figsize=(10, 6))
sns.boxplot(x='experience_level', y='salary_usd', data=df)
plt.title('Salary distribution by experience level', fontsize=14)
plt.xlabel('Experience level')
plt.ylabel('Salary (USD)')
plt.show()

# --- ANALYSIS OF PREDICTOR VARIABLE: company_size ---

plt.figure(figsize=(10, 5))

# Univariate plot
cs_count = df['company_size'].value_counts()

sns.barplot(x=cs_count.index, y=cs_count.values)
plt.title('Frequency by company size')
plt.xlabel('Company size')
plt.ylabel('Employees')

ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)

for i, v in enumerate(cs_count.values):
    plt.text(i, v + 10, str(v), ha='center', va='bottom')

plt.show()

# Bivariate plot (X vs Y)
plt.figure(figsize=(10, 6))
sns.boxplot(x='company_size', y='salary_usd', data=df)
plt.title('Salary distribution by company size', fontsize=14)
plt.xlabel('Company size')
plt.ylabel('Salary (USD)')
plt.show()

threshold = 1.5
Q1 = df['salary_usd'].quantile(0.25)
Q3 = df['salary_usd'].quantile(0.75)
IQR = Q3 - Q1

floor = Q1 - threshold * IQR
ceiling = Q3 + threshold * IQR

print(f"Interquartile Range (IQR): {IQR:,.0f}")
print(f"Tukey's fences: < {floor:,.0f} and > {ceiling:,.0f}")

Interquartile Range (IQR): 76,229
Tukey's fences: < -44,163 and > 260,752

# Identify how many points fall outside the fences
outliers = df[(df['salary_usd'] < floor) | (df['salary_usd'] > ceiling)]
print(f"Detected {len(outliers)} outliers ({len(outliers)/len(df):.2%} of the dataset).")

# Decision: filter them out
df_clean = df.drop(outliers.index).copy()

print("Treatment applied: removed records outside Tukey's fences.")

# Visualization after outlier treatment
plt.figure(figsize=(10, 4))
sns.boxplot(x=df_clean['salary_usd'])
plt.title('Distribution after removing outliers')
plt.show()

df_clean

Detected 483 outliers (3.22% of the dataset).
Treatment applied: removed records outside Tukey's fences.

# GLOBAL ANALYSIS
print("GLOBAL NORMALITY TEST (Variable: salary_usd)")
stat, p_value = stats.normaltest(df_clean['salary_usd'])

print(f"   Statistic: {stat:.4f}, p-value: {p_value:.4e}")
if p_value < 0.05:
    print("   ✅ CONCLUSION: We reject global normality (consistent with the skewness seen in EDA).")
else:
    print("   ❌ CONCLUSION: We cannot reject global normality.")

GLOBAL NORMALITY TEST (Variable: salary_usd)
   Statistic: 1287.3386, p-value: 2.8706e-280
   ✅ CONCLUSION: We reject global normality (consistent with the skewness seen in EDA).

# GROUP-BASED ANALYSIS
def check_normality_groups(df, group_col, target_col):
    # Use .cat.categories to respect our hierarchical ordering
    unique_groups = df[group_col].cat.categories

    # Plot configuration
    # Dynamically adjust the figure size based on the number of groups
    fig, axes = plt.subplots(1, len(unique_groups), figsize=(5 * len(unique_groups), 4))

    print(f"\nNORMALITY TEST BY GROUPS: {group_col}")

    for i, group in enumerate(unique_groups):
        # Filter data for the specific group
        data = df[df[group_col] == group][target_col]

        # A. Statistical test (D'Agostino)
        if len(data) >= 8:  # Minimum required for the test; otherwise it would not be valid
            stat, p_value = stats.normaltest(data)
            normality = "FAIL TO REJECT NORMALITY" if p_value > 0.05 else "NOT NORMAL"
            print(f"   -> Group '{group}' (n={len(data)}): p-value={p_value:.4e} => {normality}")
        else:
            print(f"   -> Group '{group}': Not enough data for a reliable test.")

        # B. Q-Q Plot
        ax = axes[i] if len(unique_groups) > 1 else axes
        sm.qqplot(data, line='s', ax=ax, markerfacecolor='blue', markeredgecolor='blue', alpha=0.3)
        ax.set_title(f'Q-Q Plot: {group}')

    plt.tight_layout()
    plt.show()


# Run for Experience Level
check_normality_groups(df_clean, 'experience_level', 'salary_usd')

# Run for Company Size
check_normality_groups(df_clean, 'company_size', 'salary_usd')

NORMALITY TEST BY GROUPS: experience_level
   -> Group 'EN' (n=3718): p-value=5.1094e-60 => NOT NORMAL
   -> Group 'MI' (n=3781): p-value=1.5729e-54 => NOT NORMAL
   -> Group 'SE' (n=3741): p-value=5.2348e-54 => NOT NORMAL
   -> Group 'EX' (n=3277): p-value=1.0541e-72 => NOT NORMAL

NORMALITY TEST BY GROUPS: company_size
   -> Group 'S' (n=4941): p-value=1.1475e-126 => NOT NORMAL
   -> Group 'M' (n=4855): p-value=2.7140e-108 => NOT NORMAL
   -> Group 'L' (n=4721): p-value=2.5776e-71 => NOT NORMAL

def check_homoscedasticity(df, group_col, target_col):
    # Prepare groups
    groups = [df[df[group_col] == g][target_col] for g in df[group_col].unique()]

    # Levene test (center='median' is more robust for skewed data)
    stat, p_value = stats.levene(*groups, center='median')

    print(f"--- Levene test ({group_col}) ---")
    print(f"W statistic: {stat:.4f}")
    print(f"p-value: {p_value:.4e}")

    if p_value < 0.05:
        print("✅ CONCLUSION: Variances are significantly different (heteroscedasticity).")
    else:
        print("❌ CONCLUSION: No evidence that variances differ.")

# Check for experience and company size
check_homoscedasticity(df_clean, 'experience_level', 'salary_usd')
check_homoscedasticity(df_clean, 'company_size', 'salary_usd')

--- Levene test (experience_level) ---
W statistic: 835.7112
p-value: 0.0000e+00
✅ CONCLUSION: Variances are significantly different (heteroscedasticity).
--- Levene test (company_size) ---
W statistic: 16.2998
p-value: 8.4926e-08
✅ CONCLUSION: Variances are significantly different (heteroscedasticity).

# Quantify dispersion by group to understand "who varies more" (not only who earns more)

def dispersion_table(df, group_col, target_col="salary_usd", observed=True):
    """Return a dispersion table by group (variance, std, IQR, CV)."""
    g = df.groupby(group_col, observed=observed)[target_col]
    summary = g.agg(
        n="count",
        mean="mean",
        median="median",
        std="std",
        var="var",
        q1=lambda x: x.quantile(0.25),
        q3=lambda x: x.quantile(0.75),
    )
    summary["iqr"] = summary["q3"] - summary["q1"]
    summary["cv"] = summary["std"] / summary["mean"]
    return summary.sort_values("std", ascending=False)

for col in ["experience_level", "company_size"]:
    disp = dispersion_table(df_clean, col, "salary_usd", observed=True)
    display(disp)

    g_max = disp["std"].idxmax()
    g_min = disp["std"].idxmin()
    std_ratio = disp.loc[g_max, "std"] / disp.loc[g_min, "std"]

    print(
        f"➡️ Highest dispersion (STD) in '{col}': {g_max} "
        f"(std={disp.loc[g_max,'std']:,.0f} | IQR={disp.loc[g_max,'iqr']:,.0f} | CV={disp.loc[g_max,'cv']:.2f})"
    )
    print(
        f"➡️ Lowest dispersion (STD) in '{col}': {g_min} "
        f"(std={disp.loc[g_min,'std']:,.0f} | IQR={disp.loc[g_min,'iqr']:,.0f} | CV={disp.loc[g_min,'cv']:.2f})"
    )
    print(f"📈 Dispersion ratio (std max / std min): {std_ratio:.2f}")

    # Quick visualization: standard deviation by group
    plt.figure(figsize=(8, 3))
    disp["std"].sort_values().plot(kind="barh")
    plt.title(f"Salary standard deviation by group: {col}")
    plt.xlabel("STD(salary_usd)")
    plt.ylabel(col)
    plt.show()

➡️ Highest dispersion (STD) in 'experience_level': EX (std=42,646 | IQR=64,476 | CV=0.25)
➡️ Lowest dispersion (STD) in 'experience_level': EN (std=18,558 | IQR=26,350 | CV=0.29)
📈 Dispersion ratio (std max / std min): 2.30

➡️ Highest dispersion (STD) in 'company_size': L (std=51,197 | IQR=75,071 | CV=0.43)
➡️ Lowest dispersion (STD) in 'company_size': S (std=48,243 | IQR=64,613 | CV=0.48)
📈 Dispersion ratio (std max / std min): 1.06

def test_kruskal_global(df, group_col, target_col):
    groups = [df[df[group_col] == val][target_col] for val in df[group_col].unique()]
    stat, p_value = stats.kruskal(*groups)

    print(f"Kruskal–Wallis test for variable: {group_col}")
    print(f"   H statistic: {stat:.4f} | p-value: {p_value:.4e}")

    if p_value < 0.05:
        print(f"   ✅ H0 rejected: '{group_col}' significantly affects salary.")
    else:
        print(f"   ❌ H0 not rejected: we cannot claim significant differences.")
    print("-" * 98)

# Apply to our two candidate predictors
test_kruskal_global(df_clean, 'experience_level', 'salary_usd')
test_kruskal_global(df_clean, 'company_size', 'salary_usd')

Kruskal–Wallis test for variable: experience_level
   H statistic: 9368.7891 | p-value: 0.0000e+00
   ✅ H0 rejected: 'experience_level' significantly affects salary.
--------------------------------------------------------------------------------------------------
Kruskal–Wallis test for variable: company_size
   H statistic: 449.3941 | p-value: 2.6021e-98
   ✅ H0 rejected: 'company_size' significantly affects salary.
--------------------------------------------------------------------------------------------------

# Effect size for Kruskal–Wallis: how large is the global difference?

def kruskal_effect_size(df, group_col, target_col="salary_usd", observed=True):
    groups = [g[target_col].values for _, g in df.groupby(group_col, observed=observed)]
    H, p = stats.kruskal(*groups)
    n = sum(len(x) for x in groups)
    k = len(groups)

    # Epsilon squared (ε²) is a common effect size for Kruskal–Wallis
    epsilon_sq = (H - k + 1) / (n - k)

    # Alternative: approximate eta squared (η²) for reference
    eta_sq = H / (n - 1)

    return H, p, epsilon_sq, eta_sq, n, k

for col in ["experience_level", "company_size"]:
    H, p, eps2, eta2, n, k = kruskal_effect_size(df_clean, col, "salary_usd", observed=True)
    print(f"\nEffect size for: {col}")
    print(f"H = {H:,.2f} | p = {p:.2e} | n = {n} | k = {k}")
    print(f"ε² (epsilon squared) = {eps2:.4f}")
    print(f"η² (approx.)         = {eta2:.4f}")

    # Rule of thumb (indicative): 0.01 small, 0.06 moderate, 0.14 large
    if eps2 < 0.01:
        mag = "small"
    elif eps2 < 0.06:
        mag = "moderate"
    elif eps2 < 0.14:
        mag = "medium"
    else:
        mag = "large"
    print(f"➡️ Interpretation (indicative): {mag} effect for {col}.")

Effect size for: experience_level
H = 9,368.79 | p = 0.00e+00 | n = 14517 | k = 4
ε² (epsilon squared) = 0.6453
η² (approx.)         = 0.6454
➡️ Interpretation (indicative): large effect for experience_level.

Effect size for: company_size
H = 449.39 | p = 2.60e-98 | n = 14517 | k = 3
ε² (epsilon squared) = 0.0308
η² (approx.)         = 0.0310
➡️ Interpretation (indicative): moderate effect for company_size.

def bootstrap_difference(df, group_col, control_group, test_group, target='salary_usd'):
    # Data filtering
    data_control = df[df[group_col] == control_group][target].values
    data_test = df[df[group_col] == test_group][target].values

    # Statistic function (Difference in means: Test - Control)
    # Vectorized ('np.mean(axis=-1)' is faster and more technically robust than 'data.mean()')
    def diff_means(x, y, axis=-1):
        return np.mean(y, axis=axis) - np.mean(x, axis=axis)

    # Bootstrap execution
    # We set 'random_state' for reproducibility
    res = bootstrap(
        (data_control, data_test),
        diff_means,
        confidence_level=0.95,
        random_state=24
    )

    ci = res.confidence_interval
    mean_diff = np.mean(data_test) - np.mean(data_control)

    print(f"\n📊 BOOTSTRAP: {group_col} ({test_group} vs {control_group})")
    print(f"   Observed Mean Difference: ${mean_diff:,.2f}")
    print(f"   Confidence Interval (95%): [${ci.low:,.2f}, ${ci.high:,.2f}]")

    if ci.low > 0:
        print(f"   ✅ CONCLUSION: Significant positive difference ({test_group} earns more).")
    elif ci.high < 0:
        print(f"   ✅ CONCLUSION: Significant negative difference ({control_group} earns more).")
    else:
        print("   ❌ CONCLUSION: Not significant (the interval includes 0).")


# Specific experience analysis (Senior vs Mid)
bootstrap_difference(df_clean, 'experience_level', 'MI', 'SE')

# Specific company size analysis (Medium vs Small)
# We compare M vs S because it represents the first growth jump from a startup
bootstrap_difference(df_clean, 'company_size', 'S', 'M')

📊 BOOTSTRAP: experience_level (SE vs MI)
   Observed Mean Difference: $34,232.19
   Confidence Interval (95%): [$32,857.99, $35,623.00]
   ✅ CONCLUSION: Significant positive difference (SE earns more).

📊 BOOTSTRAP: company_size (M vs S)
   Observed Mean Difference: $8,530.28
   Confidence Interval (95%): [$6,578.70, $10,509.18]
   ✅ CONCLUSION: Significant positive difference (M earns more).

# Statistical significance vs practical relevance: absolute and relative difference (%)

def pair_practical(df, group_col, group_a, group_b, target_col="salary_usd"):
    a = df.loc[df[group_col] == group_a, target_col]
    b = df.loc[df[group_col] == group_b, target_col]

    mean_a, mean_b = a.mean(), b.mean()
    med_a, med_b = a.median(), b.median()

    diff_mean = mean_a - mean_b
    pct_mean = diff_mean / mean_b * 100

    diff_med = med_a - med_b
    pct_med = diff_med / med_b * 100

    print(f"\n📌 Practical comparison: {group_col} ({group_a} vs {group_b})")
    print(f"Mean:   {group_a} = {mean_a:,.0f} | {group_b} = {mean_b:,.0f} | Δ = {diff_mean:,.0f} ({pct_mean:.1f}%)")
    print(f"Median: {group_a} = {med_a:,.0f} | {group_b} = {med_b:,.0f} | Δ = {diff_med:,.0f} ({pct_med:.1f}%)")

pair_practical(df, "experience_level", "SE", "MI", "salary_usd")
pair_practical(df, "company_size", "M", "S", "salary_usd")

📌 Practical comparison: experience_level (SE vs MI)
Mean:   SE = 122,188 | MI = 87,955 | Δ = 34,232 (38.9%)
Median: SE = 116,907 | MI = 84,641 | Δ = 32,266 (38.1%)

📌 Practical comparison: company_size (M vs S)
Mean:   M = 113,600 | S = 102,147 | Δ = 11,453 (11.2%)
Median: M = 98,122 | S = 89,648 | Δ = 8,474 (9.5%)

# TREND ANALYSIS (Spearman Correlation)
print("--- SPEARMAN CORRELATION MATRIX (Direction) ---")

# Prepare the data -> numeric encoding while preserving the hierarchy (Junior < Senior)
df_corr = df_clean.copy()
df_corr['exp_code'] = df_corr['experience_level'].cat.codes
df_corr['size_code'] = df_corr['company_size'].cat.codes

# Compute correlation (method='spearman' due to non-normality and ordinal nature; it is also more robust to outliers)
corr_matrix = df_corr[['salary_usd', 'exp_code', 'size_code']].corr(method='spearman')

# Visualization
plt.figure(figsize=(8, 6))
sns.heatmap(corr_matrix, annot=True, cmap='RdBu_r', fmt=".2f", vmin=-1, vmax=1)
plt.title("Spearman Correlation")
plt.show()

--- SPEARMAN CORRELATION MATRIX (Direction) ---

# EFFECT SIZE AND ASSOCIATION ANALYSIS (Advanced Metrics)

# First, we define the functions we will need

def compute_eta_squared(df, cat_col, num_col):
    """
    Computes effect size (Eta-squared) using the Kruskal–Wallis H statistic.
    Interpretation: % of variance in the numerical variable explained by the categorical variable.
    """
    groups = [df[df[cat_col] == val][num_col] for val in df[cat_col].unique()]
    H, _ = stats.kruskal(*groups)
    k = len(groups)     # k = number of groups
    n = len(df)         # n = total number of observations
    eta_sq = (H - k + 1) / (n - k)
    return max(0, eta_sq)

    """
    More detailed interpretation:
    - Numerator (H - k + 1)
        Adjusts the H statistic by removing the “noise” introduced by having multiple groups
    - Denominator (n - k)
        Normalizes the result so it becomes a proportion between 0 and 1
    - Result
        A value between 0 and 1, interpreted as:
            “Proportion of salary variability explained by the categorical variable”
    """


def compute_cramers_v(df, col1, col2):
    """
    Computes Cramér’s V based on the Chi-square statistic.
    Measures the association between two categorical variables (0 = none, 1 = perfect).
    """
    confusion_matrix = pd.crosstab(df[col1], df[col2])
    chi2 = stats.chi2_contingency(confusion_matrix)[0]
    n = confusion_matrix.sum().sum()

    # Normalize chi2
    phi2 = chi2 / n  # chi2 measures deviation from independence, and n = total sample size
    r, k = confusion_matrix.shape  # r = number of rows and k = number of columns

    # Bias correction (recommended for more robust results)
    with np.errstate(divide='ignore', invalid='ignore'):
        # Subtract bias to avoid inflating association
        phi2corr = max(0, phi2 - ((k - 1) * (r - 1)) / (n - 1))
        # Adjust “effective categories” to reduce error due to small cells
        rcorr = r - ((r - 1) ** 2) / (n - 1)
        kcorr = k - ((k - 1) ** 2) / (n - 1)
        # Corrected formula -> robust V between 0 and 1
        v = np.sqrt(phi2corr / min((kcorr - 1), (rcorr - 1)))

    return v

# RUN OUR ANALYSIS

print("--- EFFECT SIZE ON SALARY (Eta-squared) ---")
# Evaluate how much each predictor influences salary
for col in ['experience_level', 'company_size']:
    eta = compute_eta_squared(df_clean, col, 'salary_usd')
    print(f"Variable '{col}':")
    print(f"   -> η² = {eta:.4f} (Explains {eta*100:.1f}% of the salary variance)")

print("\n--- RELATIONSHIP BETWEEN PREDICTORS (Multicollinearity - Cramér's V) ---")
# Evaluate whether predictors provide redundant information
v_cramer = compute_cramers_v(df_clean, 'experience_level', 'company_size')
print("Association between 'experience_level' and 'company_size':")
print(f"   -> V = {v_cramer:.4f}")

# Automatic interpretation
if v_cramer < 0.3:
    print("   ✅ CONCLUSION: Low association. Predictors provide unique information (No multicollinearity risk).")
elif v_cramer < 0.7:
    print("   ⚠️ CONCLUSION: Moderate association. Watch for possible redundancy (Potential multicollinearity risk).")
else:
    print("   ❌ CONCLUSION: Very high association. Likely severe multicollinearity.")

--- EFFECT SIZE ON SALARY (Eta-squared) ---
Variable 'experience_level':
   -> η² = 0.6453 (Explains 64.5% of the salary variance)
Variable 'company_size':
   -> η² = 0.0308 (Explains 3.1% of the salary variance)

--- RELATIONSHIP BETWEEN PREDICTORS (Multicollinearity - Cramér's V) ---
Association between 'experience_level' and 'company_size':
   -> V = 0.0370
   ✅ CONCLUSION: Low association. Predictors provide unique information (No multicollinearity risk).

# VARIABLE PREPARATION (DUMMIES)

# NOTE:
# There are two common ways to do this in statsmodels:
# A) Using formulas (smf.ols): "salary_usd ~ C(experience_level) + C(company_size)"
#    More automatic because it creates dummies internally
# B) Using arrays (sm.OLS): requires creating dummies manually via pd.get_dummies
#    We'll use option (B) to keep full control over the columns and be explicit.

# Create dummy variables (one-hot encoding)
# drop_first=True removes one category as the reference (avoids perfect multicollinearity)
# dtype=int ensures 0/1 instead of True/False (avoids dtype issues in statsmodels)
X_dummies = pd.get_dummies(df_clean[['experience_level', 'company_size']], drop_first=True, dtype=int)

# Add constant (intercept) for OLS
X = sm.add_constant(X_dummies)

# Define target variable (y) and its log transform
y = df_clean['salary_usd']
y_log = np.log(df_clean['salary_usd'])  # log transform to correct heteroscedasticity

print("Variables ready. Model columns:", X.columns.tolist())
# Verify numeric types
print("\nData types in X:")
print(X.dtypes)

Variables ready. Model columns: ['const', 'experience_level_MI', 'experience_level_SE', 'experience_level_EX', 'company_size_M', 'company_size_L']

Data types in X:
const                  float64
experience_level_MI      int32
experience_level_SE      int32
experience_level_EX      int32
company_size_M           int32
company_size_L           int32
dtype: object

# PLAIN LINEAR MODEL (diagnostic baseline)

# Fit the model without log transform
model_linear = sm.OLS(y, X).fit()

print("--- LINEAR MODEL DIAGNOSTICS (NO LOG) ---")
print(f"R-squared: {model_linear.rsquared:.4f}")
print(f"Log-Likelihood: {model_linear.llf:.2f}")
print(f"AIC: {model_linear.aic:.2f}")
print(f"BIC: {model_linear.bic:.2f}")

# Residual plot to detect heteroscedasticity
residuals_linear = model_linear.resid
plt.figure(figsize=(10, 6))
plt.scatter(model_linear.fittedvalues, residuals_linear, alpha=0.5, color='gray')
plt.title("Residual diagnostics (linear model)", fontsize=14)
plt.xlabel("Predicted salary", fontsize=12)
plt.ylabel("Residuals (error)", fontsize=12)
plt.axhline(0, color='red', linestyle='--', linewidth=2)
plt.show()

--- LINEAR MODEL DIAGNOSTICS (NO LOG) ---
R-squared: 0.6522
Log-Likelihood: -170192.26
AIC: 340396.52
BIC: 340442.01

# LOG MODEL

# Fit the model using log salary
model_log = sm.OLS(y_log, X).fit()

# Show the full statistical summary
print(model_log.summary())

# ASSUMPTION DIAGNOSTICS (Gauss–Markov validation) ---
residuals = model_log.resid
fitted_vals = model_log.fittedvalues

fig, ax = plt.subplots(1, 2, figsize=(15, 6))

# -> Homoscedasticity (Residuals vs Fitted)
# We look for a shapeless "cloud". If there is no "fan" pattern, great!
sns.scatterplot(x=fitted_vals, y=residuals, ax=ax[0], alpha=0.6)
ax[0].axhline(0, color='red', linestyle='--', linewidth=2)
ax[0].set_title('Homoscedasticity', fontsize=12)
ax[0].set_xlabel('Predicted log(Salary)')
ax[0].set_ylabel('Residuals')

# -> Normality (Q-Q Plot)
# We look for blue points to follow the red line
sm.qqplot(residuals, line='45', fit=True, ax=ax[1])
ax[1].set_title('Normality (Q-Q Plot)', fontsize=12)

plt.tight_layout()
plt.show()

print("Interpretation:")
print("1. Left plot: if we see a dispersed cloud with no fan shape, we have corrected heteroscedasticity.")
print("2. Right plot: if the points align with the red line, residuals are approximately normally distributed.")

                            OLS Regression Results                            
==============================================================================
Dep. Variable:             salary_usd   R-squared:                       0.679
Model:                            OLS   Adj. R-squared:                  0.678
Method:                 Least Squares   F-statistic:                     6127.
Date:                Wed, 14 Jan 2026   Prob (F-statistic):               0.00
Time:                        21:35:05   Log-Likelihood:                -1170.5
No. Observations:               14517   AIC:                             2353.
Df Residuals:                   14511   BIC:                             2399.
Df Model:                           5                                         
Covariance Type:            nonrobust                                         
=======================================================================================
                          coef    std err          t      P>|t|      [0.025      0.975]
---------------------------------------------------------------------------------------
const                  10.9016      0.005   2071.176      0.000      10.891      10.912
experience_level_MI     0.3288      0.006     54.258      0.000       0.317       0.341
experience_level_SE     0.6576      0.006    108.225      0.000       0.646       0.669
experience_level_EX     1.0173      0.006    161.743      0.000       1.005       1.030
company_size_M          0.0998      0.005     18.818      0.000       0.089       0.110
company_size_L          0.2345      0.005     43.853      0.000       0.224       0.245
==============================================================================
Omnibus:                     1233.878   Durbin-Watson:                   1.986
Prob(Omnibus):                  0.000   Jarque-Bera (JB):              396.771
Skew:                           0.067   Prob(JB):                     6.95e-87
Kurtosis:                       2.201   Cond. No.                         5.29
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Interpretation:
1. Left plot: if we see a dispersed cloud with no fan shape, we have corrected heteroscedasticity.
2. Right plot: if the points align with the red line, residuals are approximately normally distributed.

# ADDITIONAL RESIDUAL PLOTS (Histogram and Boxplot)
fig, ax = plt.subplots(1, 2, figsize=(15, 6))

# Histogram
sns.histplot(residuals, bins=30, kde=True, ax=ax[0])
ax[0].set_title("Residuals histogram", fontsize=12)
ax[0].set_xlabel("Error (residuals)")

# Boxplot (to see outliers)
sns.boxplot(x=residuals, ax=ax[1])
ax[1].set_title("Residuals boxplot", fontsize=12)
ax[1].set_xlabel("Error distribution")

plt.show()

# Shape metrics
print(f"Skewness: {residuals.skew():.4f}")
print(f"Kurtosis: {residuals.kurtosis():.4f}")

Skewness: 0.0673
Kurtosis: -0.7985

# A) PREDICTION FOR A HYPOTHETICAL CASE

# Let's imagine a new case: a 'Senior' (SE) in a 'Large' (L) company
new_employee = pd.DataFrame(columns=X.columns)
new_employee.loc[0] = 0  # initialize all predictors to 0
new_employee['const'] = 1

# Activate the corresponding dummy variables
if 'experience_level_SE' in X.columns:
    new_employee['experience_level_SE'] = 1
if 'company_size_L' in X.columns:
    new_employee['company_size_L'] = 1

# Prediction and inverse log-transform (np.exp)
pred_log = model_log.predict(new_employee)[0]
pred_usd = np.exp(pred_log)

# Analytical confidence interval (95%)
pred_summary = model_log.get_prediction(new_employee).summary_frame(alpha=0.05)
ci_lower = np.exp(pred_summary['obs_ci_lower'][0])
ci_upper = np.exp(pred_summary['obs_ci_upper'][0])

print("--- CASE PREDICTION: Senior in a Large Company ---")
print(f"Estimated salary: ${pred_usd:,.2f}")
print(f"Interval (95%): [${ci_lower:,.2f} - ${ci_upper:,.2f}]")


# B) ROBUST VALIDATION (BOOTSTRAPPING)

rmse_list = []
n_iter = 2000

# Set a seed for reproducibility
np.random.seed(24)

print(f"\nRunning bootstrapping ({n_iter} iterations)...")

for _ in tqdm(range(n_iter)):
    # Split into train/test (randomness is controlled by the fixed seed)
    X_train, X_test, y_train, y_test = train_test_split(X, y_log, test_size=0.2)  # note: no random_state here

    # Train and predict
    m_boot = sm.OLS(y_train, X_train).fit()  # model is fitted here
    preds_log = m_boot.predict(X_test)

    # Convert back to real USD to compute the true error
    rmse = np.sqrt(mean_squared_error(np.exp(y_test), np.exp(preds_log)))
    rmse_list.append(rmse)

# Compute mean and percentiles
rmse_mean = np.mean(rmse_list)
rmse_lower = np.percentile(rmse_list, 2.5)
rmse_upper = np.percentile(rmse_list, 97.5)

print(f"Global mean RMSE: ${rmse_mean:,.2f}")
print(f"Error range (95%): [${rmse_lower:,.0f} - ${rmse_upper:,.0f}]")

--- CASE PREDICTION: Senior in a Large Company ---
Estimated salary: $132,412.64
Interval (95%): [$79,169.73 - $221,462.28]

Running bootstrapping (2000 iterations)...

100%|██████████| 2000/2000 [00:18<00:00, 108.63it/s]

Global mean RMSE: $30,146.95
Error range (95%): [$29,360 - $30,976]

# Prepare a sample and create numeric codes on the fly (also set a seed)
# We use sample(500) to visualize density without overcrowding the plot
df_viz = df_clean.sample(500, random_state=24).copy()

# Since we defined the logical ordering (Entry -> Executive) at the beginning of the analysis,
# we can use .cat.codes to automatically respect that ordering
df_viz['exp_code_viz'] = df_viz['experience_level'].cat.codes
df_viz['size_code_viz'] = df_viz['company_size'].cat.codes

# Generate predictions for this sample
# Recreate the dummy variables for the model
X_viz = pd.get_dummies(df_viz[['experience_level', 'company_size']], drop_first=True, dtype=int)
X_viz = sm.add_constant(X_viz)

# Fill missing columns with 0 (in case the sample is missing a category)
for col in X.columns:
    if col not in X_viz.columns:
        X_viz[col] = 0
X_viz = X_viz[X.columns]  # reorder columns to match the original model

# Predict
pred_log_viz = model_log.predict(X_viz)
df_viz['pred_usd'] = np.exp(pred_log_viz)

# Build the INTERACTIVE 3D PLOT
fig = go.Figure()

# A) Actual points (Blue)
fig.add_trace(go.Scatter3d(
    x=df_viz['exp_code_viz'],
    y=df_viz['size_code_viz'],
    z=df_viz['salary_usd'],
    mode='markers',
    marker=dict(size=5, color='blue', opacity=0.7),
    name="Actual Data"
))

# B) Predicted points (Orange)
fig.add_trace(go.Scatter3d(
    x=df_viz['exp_code_viz'],
    y=df_viz['size_code_viz'],
    z=df_viz['pred_usd'],
    mode='markers',
    marker=dict(size=5, color='orange', opacity=0.9, symbol='diamond'),
    name="Model Prediction"
))

# C) Residual lines (Red)
for x_real, y_real, z_real, z_pred in zip(
    df_viz['exp_code_viz'],
    df_viz['size_code_viz'],
    df_viz['salary_usd'],
    df_viz['pred_usd']
):
    fig.add_trace(go.Scatter3d(
        x=[x_real, x_real],
        y=[y_real, y_real],
        z=[z_real, z_pred],
        mode='lines',
        line=dict(color='red', width=2),
        showlegend=False  # hide from the legend
    ))

# Final layout configuration with axis text labels
fig.update_layout(
    title="Prediction Errors: Actual vs Model (Sample of 500 cases)",
    scene=dict(
        xaxis=dict(title="Experience Level", tickvals=[0, 1, 2, 3], ticktext=["Entry", "Mid", "Senior", "Executive"]),
        yaxis=dict(title="Company Size", tickvals=[0, 1, 2], ticktext=["Small", "Medium", "Large"]),
        zaxis_title="Salary ($)"
    ),
    width=900,
    height=700
)

fig.show()

# Convert the bootstrapped error list to a DataFrame
errors_df = pd.DataFrame({"rmse": rmse_list})

# RMSE histogram
plt.figure(figsize=(10, 6))
errors_df['rmse'].hist(bins=30, edgecolor='black', alpha=0.7)
plt.title("Distribution of mean error (RMSE) across 2000 simulations")
plt.xlabel("Mean error in dollars ($)")
plt.ylabel("Frequency")
plt.axvline(errors_df['rmse'].mean(), color='red', linestyle='dashed', linewidth=2, label='Mean')
plt.legend()
plt.show()

# Distribution statistics
print(f"Skewness of error: {errors_df['rmse'].skew():.4f}")
print(f"Kurtosis of error: {errors_df['rmse'].kurtosis():.4f}")

# Confidence interval
mean_error = errors_df['rmse'].mean()
std_error = errors_df['rmse'].std()

lower = np.percentile(errors_df['rmse'], 2.5)
upper = np.percentile(errors_df['rmse'], 97.5)

print("\n--- FINAL ACCURACY CONCLUSION (Percentile method) ---")
print(f"The model's mean RMSE is ${mean_error:,.0f}.")
print("With 95% confidence, the prediction error will lie within:")
print(f"[{lower:,.0f}, {upper:,.0f}]")
print(f"This shows high model stability, with a standard deviation of only ${std_error:,.0f} across simulations.")

Skewness of error: 0.0438
Kurtosis of error: 0.0910

--- FINAL ACCURACY CONCLUSION (Percentile method) ---
The model's mean RMSE is $30,147.
With 95% confidence, the prediction error will lie within:
[29,360, 30,976]
This shows high model stability, with a standard deviation of only $412 across simulations.

Field Name	Description	Data Type	Format / Example	Possible Values / Range	Nulls Allowed	Notes	Semantic Type
`job_id`	Unique identifier for each job posting	Categorical (string)	'AI00001'	N/A	No	Primary key, no duplicates.	Categorical
`job_title`	Job title or role advertised	Categorical (string)	'AI Research Scientist'	20	No	Includes all possible roles. 20 unique values.	Categorical
`salary_usd`	Annual salary converted to USD	Integer (int)	90,376	≥ 0	No	Salary expressed in US dollars.	Ratio
`salary_currency`	Currency in which the salary is offered	Categorical (string)	'USD'	'USD', 'EUR', 'GBP'	No	USD = US Dollar, EUR = Euro, GBP = Pound Sterling.	Categorical
`experience_level`	Required level of professional experience	Categorical (string)	'SE'	'EN', 'MI', 'SE', 'EX'	No	EN = Entry, MI = Mid, SE = Senior, EX = Executive.	Categorical / Ordinal
`employment_type`	Type of employment contract	Categorical (string)	'CT'	'FT', 'PT', 'CT', 'FL'	No	CT = Contract, FT = Full-time, PT = Part-time, FL = Freelance.	Categorical
`company_location`	Country where the company is based	Categorical (string)	'China'	20	No	20 unique country values.	Categorical
`company_size`	Size of the company	Categorical (string)	'S'	3	No	S = Small, M = Medium, L = Large.	Categorical / Ordinal
`employee_residence`	Country of employee residence	Categorical (string)	'China'	20	No	20 unique country values.	Categorical
`remote_ratio`	Percentage of remote work allowed	Integer (int)	50	{0, 50, 100}	No	0 = On-site, 50 = Hybrid, 100 = Fully remote.	Ratio
`required_skills`	Key technical skills required (comma-separated)	Categorical (string)	'Python'	13,663	No	Skills are comma-separated. 13,663 unique values.	Categorical
`education_required`	Minimum education level required	Categorical (string)	'Associate'	'Associate', 'Bachelor', 'Master', 'PhD'	No	4 unique education levels.	Categorical / Ordinal
`years_experience`	Minimum years of experience required	Integer (int)	4	0 – 20	No	Discrete numeric values.	Ratio
`industry`	Industry sector of the job	Categorical (string)	'Media'	15	No	15 unique industry sectors.	Categorical
`posting_date`	Date when the job was posted	Date	'2024-10-18'	N/A	No	ISO format: YYYY-MM-DD.	Interval
`application_deadline`	Final application deadline	Date	'2024-11-07'	N/A	No	ISO format: YYYY-MM-DD.	Interval
`job_description_length`	Length of the job description (characters)	Integer (int)	1,076	N/A	No	Number of characters in the description text.	Ratio
`benefits_score`	Numeric score reflecting benefits quality	Float	9.4	0 – 10	No	Higher values indicate better benefits.	Ratio
`company_name`	Name of the hiring company	Categorical (string)	'Smart Analytics'	16	No	16 unique company names.	Categorical

	count	mean	min	25%	50%	75%	95%	99%	max	std
salary_usd	15000.0	115348.965133	32519.0	70179.75	99705.0	146408.5	237987.95	307775.69	399095.0	60260.940438
remote_ratio	15000.0	49.483333	0.0	0.0	50.0	100.0	100.0	100.0	100.0	40.812712
years_experience	15000.0	6.2532	0.0	2.0	5.0	10.0	17.0	19.0	19.0	5.545768
posting_date	15000	2024-08-29 08:48:51.840000	2024-01-01 00:00:00	2024-04-29 00:00:00	2024-08-28 00:00:00	2024-12-29 00:00:00	2025-04-06 00:00:00	2025-04-26 00:00:00	2025-04-30 00:00:00	NaN
application_deadline	15000	2024-10-11 21:55:23.520000	2024-01-16 00:00:00	2024-06-13 00:00:00	2024-10-12 00:00:00	2025-02-10 00:00:00	2025-05-21 00:00:00	2025-06-19 00:00:00	2025-07-11 00:00:00	NaN
job_description_length	15000.0	1503.314733	500.0	1003.75	1512.0	2000.0	2402.0	2479.0	2499.0	576.127083
benefits_score	15000.0	7.504273	5.0	6.2	7.5	8.8	9.8	9.9	10.0	1.45087

	n	mean	median	std	var	q1	q3	iqr	cv
experience_level
EX	3277	171357.764724	167899.0	42645.767697	1.818662e+09	138802.0	203278.0	64476.0	0.248870
SE	3741	122187.657845	116907.0	35261.617316	1.243382e+09	94173.0	145315.0	51142.0	0.288586
MI	3781	87955.471833	84641.0	25543.109400	6.524504e+08	67621.0	104483.0	36862.0	0.290410
EN	3718	63133.377084	60373.5	18558.096806	3.444030e+08	48516.0	74866.5	26350.5	0.293951

	n	mean	median	std	var	q1	q3	iqr	cv
company_size
L	4721	120178.613429	111160.0	51196.951526	2.621128e+09	80483.0	155554.0	75071.0	0.426007
M	4855	108280.711843	96759.0	50461.054640	2.546318e+09	70841.5	139181.0	68339.5	0.466021
S	4941	99750.428658	88895.0	48242.994196	2.327386e+09	63454.0	128067.0	64613.0	0.483637

	job_id	job_title	salary_usd	salary_currency	experience_level	employment_type	company_location	company_size	employee_residence	remote_ratio	required_skills	education_required	years_experience	industry	posting_date	application_deadline	job_description_length	benefits_score	company_name
0	AI00001	AI Research Scientist	90376	USD	SE	CT	China	M	China	50	Tableau, PyTorch, Kubernetes, Linux, NLP	Bachelor	9	Automotive	2024-10-18	2024-11-07	1076	5.9	Smart Analytics
1	AI00002	AI Software Engineer	61895	USD	EN	CT	Canada	M	Ireland	100	Deep Learning, AWS, Mathematics, Python, Docker	Master	1	Media	2024-11-20	2025-01-11	1268	5.2	TechCorp Inc
2	AI00003	AI Specialist	152626	USD	MI	FL	Switzerland	L	South Korea	0	Kubernetes, Deep Learning, Java, Hadoop, NLP	Associate	2	Education	2025-03-18	2025-04-07	1974	9.4	Autonomous Tech
3	AI00004	NLP Engineer	80215	USD	SE	FL	India	M	India	50	Scala, SQL, Linux, Python	PhD	7	Consulting	2024-12-23	2025-02-24	1345	8.6	Future Systems
4	AI00005	AI Consultant	54624	EUR	EN	PT	France	S	Singapore	100	MLOps, Java, Tableau, Python	Master	0	Media	2025-04-15	2025-06-23	1989	6.6	Advanced Robotics

	Column	Non-Null	Python Type	Num. Unique Values	Sample Values
0	job_id	15000	object	15000	[AI00001, AI00002, AI00003, AI00004, AI00005]
1	job_title	15000	object	20	[AI Research Scientist, AI Software Engineer, ...
2	salary_usd	15000	int64	14315	[90376, 61895, 152626, 80215, 54624]
3	salary_currency	15000	object	3	[USD, EUR, GBP]
4	experience_level	15000	object	4	[SE, EN, MI, EX]
5	employment_type	15000	object	4	[CT, FL, PT, FT]
6	company_location	15000	object	20	[China, Canada, Switzerland, India, France]
7	company_size	15000	object	3	[M, L, S]
8	employee_residence	15000	object	20	[China, Ireland, South Korea, India, Singapore]
9	remote_ratio	15000	int64	3	[50, 100, 0]
10	required_skills	15000	object	13663	[Tableau, PyTorch, Kubernetes, Linux, NLP, Dee...
11	education_required	15000	object	4	[Bachelor, Master, Associate, PhD]
12	years_experience	15000	int64	20	[9, 1, 2, 7, 0]
13	industry	15000	object	15	[Automotive, Media, Education, Consulting, Hea...
14	posting_date	15000	datetime64[ns]	486	[2024-10-18 00:00:00, 2024-11-20 00:00:00, 202...
15	application_deadline	15000	datetime64[ns]	543	[2024-11-07 00:00:00, 2025-01-11 00:00:00, 202...
16	job_description_length	15000	int64	2000	[1076, 1268, 1974, 1345, 1989]
17	benefits_score	15000	float64	51	[5.9, 5.2, 9.4, 8.6, 6.6]
18	company_name	15000	object	16	[Smart Analytics, TechCorp Inc, Autonomous Tec...

Statistical Analysis of AI Salaries (2025) 🤖

Table of Contents¶

1. Data Loading and Setup¶

2. Exploratory Data Analysis (EDA)¶

2.1 Initial Review¶

2.2 Visualizing Variable Distributions¶

2.2.1 Target Variable¶

2.2.2 Predictor Variables¶

3. Outliers¶

3.1 Identification¶

3.2 Treatment¶

4. Statistical Inference¶

4.1 Normality Assessment¶

Interpretation: why “non-normal” and what does it imply?¶

4.2 Homoscedasticity Assessment (Equal Variances)¶

Interpretation: what does it mean if variances are different?¶

Interpreting the dispersion table: who varies more and why does it matter?¶

4.3 Overall Analysis: Kruskal–Wallis¶

Assumptions check (parametric requirements)¶

Why Kruskal–Wallis?¶

Hypotheses (Kruskal–Wallis)¶

Effect size: how large is the difference?¶

4.4 Specific Analysis (Post-hoc): Quantifying Differences (Bootstrapping)¶

5. Correlation and Association Analysis¶

6. Predictive Modeling: Multiple Linear Regression¶

6.1 Variable Preparation: Dummy Encoding¶

6.2 Model Selection: Linear vs. Logarithmic¶

Heteroscedasticity Correction: Logarithmic Transformation¶

6.3 Detailed Interpretation of Statistical Results¶

1. Goodness of Fit and Information Criteria¶

2. Economic Interpretation of the Coefficients ($\beta$)¶

3. Assumption Validation (Residual Diagnostics)¶

Modeling Conclusions¶

6.4 Prediction and Confidence Intervals¶

6.5 3D Visualization (Error Visualization)¶

6.6 Model Uncertainty Analysis¶

7. Final Conclusions¶

1. Hierarchy of Salary Drivers¶

2. Market Behavior¶

3. Reliability of the Predictive Model¶

4. Executive Summary¶